#include "./dynamics.h"

void  mds_dynamics_full(mds_frame         *frame,
						mds_potential_proc Vproc, 
						void              *Vargs,
						void              *param)
{
#if defined(__INTEL_COMPILER)
#	pragma warning ( disable : 869 )
#endif
	assert( frame );
	assert( frame->sets );
	assert( frame->num_sets >= 3 );
	assert( Vproc );
	assert( frame->pbc_for_r );
	assert( frame->pbc_for_f );
	{
		const mds_vector *   position  = frame->sets[MDS_POSITION];
		mds_vector       *   force     = frame->sets[MDS_ACCELERATION];
		const size_t         n         = frame->num_particles;
		const mds_vector    *box       = & frame->box;
		mds_pbc_proc         pbc_for_f = frame->pbc_for_f;
		mds_vector           f_ij;
		mds_vector           r_ij;
		Real                 V         = 0;
		Real                 W         = 0;
		register size_t      i,j;
		
		/** initialize forces **/
		mds_frame_no_forces(frame);
		
		/******************/
		/** outer loop  **/
		/****************/
		for( i=0; i < n; ++i )
		{
			const mds_vector *r_i = &position[i];
			mds_vector       *f_i = &force[i];
			/****************/
			/** inner loop */
			/**************/
			for( j=i+1; j < n; ++j )
			{
				const mds_vector *r_j = &position[j];
				mds_vector       *f_j = &force[j];
				
				/** compute pair separation */
				mds_vector_vec( &r_ij, r_i, r_j );
				
				/** apply PBC for pair separation */
				pbc_for_f( &r_ij, box );
				
				/** update V and compute f_ij */
				V += Vproc( &r_ij, &f_ij, Vargs );
				/*
				 fprintf( stderr, "r=[%g %g %g] => V=%g\n", r_ij.x, r_ij.y, r_ij.z, Vproc( &r_ij, &f_ij, Vargs) );
				 */
				
				/** dispatch forces */
				mds_vector_sub( f_i, &f_ij );
				mds_vector_add( f_j, &f_ij );
				
				/** update Virial */
				W += mds_vector_dot( &r_ij, &f_ij );
			}
		}
		
		/** Update frame info */
		frame->V =  V;
		frame->W =  W/MDS_F(3.0);
		
	}
	
}

void mds_forces_to_accelerations( mds_frame *frame )
{
	assert( frame );
	assert( frame->sets );
	assert( frame->num_sets >= 3 );
	{
		const mds_vector    *position  = frame->sets[MDS_POSITION];
		mds_vector          *force     = frame->sets[MDS_ACCELERATION];
		const size_t         n         = frame->num_particles;
		register size_t      i;
		
		/** compute accelerations */
		for( i=0; i < n; ++i )
		{
			const mds_vector *r_i = & position[i];
			mds_vector       *f_i = & force[i];
			
			mds_vector_mul(f_i, r_i->info.inv_mass);
			
		}
		
	}
}

void mds_dynamics_call( mds_dynamics *dynamics, mds_frame *frame, mds_potential_proc Vproc, void *Vargs )
{
	assert( dynamics );
	assert( frame    );
	assert( Vproc    );
	assert( dynamics->gather );
	
	/** compute forces **/
	dynamics->gather( frame, Vproc, Vargs, dynamics->params );
	
	/** compute accelerations **/
	mds_forces_to_accelerations( frame );
}

void mds_dynamics_full_simple( mds_frame         *frame, 
							  mds_potential_proc Vproc, 
							  void              *Vargs,
							  void              *params)
{
#if defined(__INTEL_COMPILER)
#	pragma warning ( disable : 869 )
#endif
	assert( frame );
	assert( frame->sets );
	assert( frame->num_sets >= 3 );
	assert( Vproc );
	assert( frame->pbc_for_r );
	assert( frame->pbc_for_f );
	{
		const mds_vector *   position  = frame->sets[MDS_POSITION];
		mds_vector       *   force     = frame->sets[MDS_ACCELERATION];
		const size_t         n         = frame->num_particles;
		const mds_vector    *box       = & frame->box;
		mds_pbc_proc         pbc_for_f = frame->pbc_for_f;
		mds_vector           f_ij;
		mds_vector           r_ij;
		register size_t      i,j;
		
		/** initialize forces **/
		mds_frame_no_forces(frame);
		
		/******************/
		/** outer loop  **/
		/****************/
		for( i=0; i < n; ++i )
		{
			const mds_vector *r_i = &position[i];
			mds_vector       *f_i = &force[i];
			/****************/
			/** inner loop */
			/**************/
			for( j=i+1; j < n; ++j )
			{
				const mds_vector *r_j = &position[j];
				mds_vector       *f_j = &force[j];
				
				/** compute pair separation */
				mds_vector_vec( &r_ij, r_i, r_j );
				
				/** apply PBC for pair separation */
				pbc_for_f( &r_ij, box );
				
				/** update V and compute f_ij */
				Vproc( &r_ij, &f_ij, Vargs );
				
				/** dispatch forces */
				mds_vector_sub( f_i, &f_ij );
				mds_vector_add( f_j, &f_ij );
				
			}
		}
	}
}
